# 参考文献：https://blog.csdn.net/qq_44648285/article/details/143691462
import numpy as np
import matplotlib.pyplot as plt

def mvdr_beamforming(X, a_theta0):
    """
    MVDR波束成形
    X: 接收信号矩阵 (N x K)
    a_theta0: 期望方向的阵列响应向量 (N x 1)
    """
    R = np.dot(X, X.conj().T) / X.shape[1]  # 计算协方差矩阵
    R_inv = np.linalg.inv(R)
    numerator = np.dot(R_inv, a_theta0)
    denominator = np.dot(a_theta0.conj().T, numerator)
    w_mvdr = numerator / denominator
    return w_mvdr

# 示例参数
N = 8          # 天线数
K = 1000       # 快照数
theta0 = 30    # 期望信号方向（度）
theta1 = 60    # 干扰信号方向（度）
signal_power = 1
interference_power = 0.5
noise_power = 0.1

# 生成阵列响应向量
def steering_vector(theta, N, d=0.5, wavelength=1):
    angles = np.deg2rad(theta)
    return np.exp(-1j * 2 * np.pi * d * np.arange(N) * np.sin(angles) / wavelength)

a_theta0 = steering_vector(theta0, N) # 期望信号
a_theta1 = steering_vector(theta1, N) # 干扰信号

# 生成信号
s = np.sqrt(signal_power) * (np.random.randn(K) + 1j * np.random.randn(K))
i = np.sqrt(interference_power) * (np.random.randn(K) + 1j * np.random.randn(K))
n = np.sqrt(noise_power) * (np.random.randn(N, K) + 1j * np.random.randn(N, K))

# 接收信号矩阵
X = np.outer(a_theta0, s) + np.outer(a_theta1, i) + n

# 计算MVDR权重
w = mvdr_beamforming(X, a_theta0)

# 输出信号
y = np.dot(w.conj().T, X)

# 绘制信号功率谱
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
plt.figure(figsize=(10, 6))
angles = np.linspace(-90, 90, 181)
responses = []
for theta in angles:
    a_theta = steering_vector(theta, N)
    response = np.abs(np.dot(w.conj().T, a_theta))**2
    responses.append(response)
plt.plot(angles, 10 * np.log10(responses / np.max(responses)), label='MVDR Beam Pattern')
plt.title('MVDR波束图')
plt.xlabel('方向角度（度）')
plt.ylabel('归一化增益（dB）')
plt.grid(True)
plt.legend()
plt.show()

# 输出权重向量
print("MVDR权重向量:")
print(w)
